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We present an analytical model of a single natural guide star astronomical adaptive optics 
system, in closed loop mode. The model is used to simulate the long exposure system point 
spread function, using the spatial frequency (or Fourier) approach, and complement an 
initial open loop model. Applications range from system design, science case analysis and 
^ \ AO data reduction. All the classical phase errors have been included: deformable mirror 
J> ! fitting error, wavefront sensor spatial aliasing, wavefront sensor noise, and the correlated 
anisoplanatic and servo-lag error. The model includes the deformable mirror spatial 
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transfer function, and the actuator array geometry can be different from the wavefront 
sensor lenslet array geometry. We also include the dispersion between the sensing and the 
correction wavelengths. Illustrative examples are given at the end of the paper. 
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j§; 1 INTRODUCTION 

Astronomical adaptive optics (AO) are complex opto-electro-mechanical systems, de- 

signed to correct random aberrations generated by optical turbulence in earth's atmo- 
00 . 

iy) . sphere, and improve the performance of astronomical telescopes - see Roddier et al. [1] 
and figure [T] for an illustration. Modeling the performance of such a system - for science 
O ' cases analysis, AO data reduction and system design - requires sophisticated simulations 
tools. The most intuitive approach to construct an AO simulation tool is to break-down 
the system in its fundamental components, build a physical model for each of these compo- 
nents, and link all of them following a block-diagram architecture. The turbulent phase is 
then propagated through the model, mimicking a real system, up to the system's output, 
the focal plane, where the system's point spread function (PSF) is measured. The capac- 
ity of such end-to-end models to produce accurate performance predictions is in principle 
limited only by our capability to accurately code the behavior of each sub components 
and the system's input disturbance (optical turbulence, noise, other) - see for instance 
Carbillet et al. [2] and Le Louarn et al. [3]. 

End-to-end models have a limitation, though. Indeed, as the input of the system (optical 
turbulence) is of stochastic nature, the model needs to be run on a large number of 
instantaneous turbulent optical waves, typically more than a thousand to converge to a 
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Figure 1: Top: basic elements of an AO closed loop system. The incoming turbulent 
wavefront is transmitted by the telescope optics to the entrance of the AO system. The 
deformable mirror surface shape (DM) is set by the control computer to compensate for 
the wavefront error. A fraction of the corrected light beam is transmitted by the beam- 
splitter (BS) to the wavefront sensor (WFS) where the residual wavefront is measured. 
The computer reconstruct the residual wavefront from the WFS measurement and com- 
putes/updates the DM surface to keep the wavefront residual as small as possible. The 
other part of the corrected beam is sent toward the science instrument. Bottom: equiv- 
alent block diagram of the AO closed loop system, indicating the main elements of the 
loop. R indicates the wavefront reconstruction operation, C is the control algorithm (an 
integrator in this paper). 
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pseudo-long exposure PSF from which performance metrics can be estimated. Such a 
procedure takes a lot of time: depending on order of the AO system, which defines the 
size of the command matrix, and the level of sophistication of the end-to-end model, 
getting a PSF without too many residual speckles can take several hours or even days. As 
a consequence, performance analysis using end-to-end tools is generally limited to a few 
well selected representative cases, and extensive studies of the system parameter space is 
rarely undertaken. Also, because all the error sources are naturally merged in an end- 
to-end model, it is not easy to disentangle the impact of the individual sources of error 
on the overall performance, unless one has a deep understanding of how the end-to-end 
model was built. End-to-end models are therefore not a good choice for AO engineers for 
the broad analysis of a system performance, nor for astronomers interested in exploring 
the capability of a given AO system (existent or planned) for their science programme. 

To suppress the limitation of end-to-end models, Rigaut et al. jl] proposed a totally 
different method, that we refer to here as the synthetic approach. Instead of letting the 
optical wave propagate through the system components and observe what comes at the 
output, we built a model for the system's output itself. The synthetic approach is based 
on an understanding of the system's behavior, and its accuracy is only limited by this 
understanding. A priori knowledge is critical for the synthetic approach, while this is 
not needed with the end-to-end approach. We can also said that the synthetic approach 
models the behavior of the system, while the end-to-end approach models the structure 
of the system. 

The starting/central point in the synthetic approach is the construction of an analyti- 
cal model for the long exposure (or average) AO-corrected phase spatial power spectrum 
density (s-PSD). AO correction is actually seen as a spatial filter applied on the turbu- 
lent phase s-PSD. This approach is therefore also called the Fourier method or spatial 
frequency method in the AO literature. From this residual phase s-PSD it is shown in 
this paper how the long exposure PSF can be computed, in a few steps. Getting the 
long exposure PSF is therefore very fast and exploring in detail the AO system parameter 
space becomes possible. Also, wavefront error budgets are easily build, because a s-PSD 
is attributed to each error source, from which we can get the wavefront error variance, by 
numerical integration in the spatial frequency domain. Finally, by nature, there are no 
residual speckles in a synthetic PSF: performance metrics (Strehl ratio, PSF width, inte- 
grated energy) are therefore not " noisy" , and science case analysis are greatly facilitated. 

This being said, the synthetic approach has its limits: non linear effects cannot be mod- 
eled, neither transient temporal behaviors. Also, the fundamental assumption is that the 
corrected phase is stationary within the pupil (which is not true near the edge of the pupil 
and for low order aberrations) and this can produce pessimistic performance predictions. 
For these reasons, end-to-end and synthetic models are to be seen as complementary 
rather than competitive methods: broad exploration of the system's performance is the 
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domain of synthetic models, while detailed analysis of specific aspects of the system def- 
initely requires end-to-end models. Actually, inclusion of our synthetic model within an 
end-to-end Monte-Carlo code has already been tried successfully by Carbillet et al. [5]. 

In an earlier work [B], we complemented Rigaut et al.'s work, explaining the foundations 
of the method in great details, and including PSF modeling in two dimensions. This 
initial model was open loop, where the wavefront sensor (WFS) measures the wavefront 
aberration before the deformable mirror (DM) correction, while the vast majority of cur- 
rent systems operate in a closed loop mode, i.e. the WFS measures the residual wavefront 
error after compensation by the DM - see figure [TJ 

Closed loop systems control the wavefront error and are therefore relatively insensitive 
to external disturbances (like noise) and internal variations of system parameters. Open 
loop systems on the contrary are very sensitive to modeling errors: it is critical that 
the system's components behave the way they are supposed too, as the quality of the 
correction is not controlled. On the other hand, feedback of the wavefront error in a 
closed loop system can generate diverging instabilities, if the error is overcompensated, or 
if time delay in the loop is too large. Open loop systems do not have this stability issue. 
See Ogata [7] for an introduction on control systems. 

In bright guide star conditions, the AO system can be run at a high loop rate: the servo 
lag error (due to the time lag between the measurement and the actual correction) is low, 
and the noise level is negligible. In this case there are no differences between open and 
closed loop performance. For a dim guide star, the WFS exposure time is increased to 
gather more photons and keep the signal/noise ratio of the wavefront error measurement 
at an acceptable level: servo-lag error increases, and because of differences in open and 
closed loop transfer functions, the system performance significantly changes between the 
two modes. Because of this different behavior, open loop models cannot be used to predict 
close loop performance in high noise regime. As predicting the limiting magnitude of a 
system is of great importance, in particular for science cases studies, we have developed 
further our initial model and included closed loop modeling. 

We also took this opportunity to introduce a DM spatial frequency transfer function, 
allowing the analysis of different influence function structures and actuators grid architec- 
ture, and the dispersion error, i.e the error induced by the air's refractive index dispersion, 
which makes that the wavefront measured at a given wavelength is slightly different from 
the wavefront corrected at any other wavelength, generating a non negligible error for 
tight error budget AO systems. Also, a few conceptual errors that appeared in our initial 
open loop paper are corrected. 

Our closed loop model is developed for a single natural guide star (NGS) Shack-Hartmann 
WFS based system (SH-WFS), and a least square error (LSE) wavefront reconstruction 
algorithm. This case covers the vast majority of current NGS-based AO systems design, 
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and in any case our model can be easily adapted to other schemes. Examples of usage of 
the synthetic method are given at the end of this paper, and we show how this method 
can be used to optimally dimension a system. 

We finally note that the synthetic approach has been used and developed by other authors 
as well, increasing the diversity of views and understandings of the limits and strengths 
of this method. We must mention first the work of Ellerbroek [8] who basically developed 
the same method using a more general albeit potentially less detailed approach; the work 
of Rigaut [9], Tokovinin [10] and Jolissaint et al. [11] for a ground-layer AO mode (but not 
in closed loop and without wavefront sensor noise model); and more recently Neichel et 
al. [12] who introduced multiple guide star tomographic reconstruction, a very useful ex- 
tension of the method. What our model brings to these recent developments is essentially 
the closed loop mode, and some useful sophistications like the DM transfer function. To 
finish, it is fair to mention that R. Conan and Ch. Verinaud (private communication) both 
independently developed a synthetic closed loop model, yet unpublished, using another 
but equivalent approach than the one presented here, namely the equivalence between the 
spatial and temporal frequency through the Taylor hypothesis (see text). 

2 FOUNDATIONS OF THE SPATIAL FREQUENCY 
METHOD: A SUMMARY 

A detailed description of the foundations of the method is given in Jolissaint et al. [B] . A 
summary is given here for convenience. 

The method is based on the relationship between the phase spatial frequency power 
spectrum and the phase structure function (SF) in one side, and the SF and the average 
(long exposure) AO system optical transfer function (OTF) on the otheJl Let us review 
the different steps from the phase s-PSD to the long exposure PSF. 

Analytical expressions for the s-PSD will be developed later. The starting point is the 
OTF of the whole system made of (1) the column of air above the telescope, seen as an 
optically transparent medium with a turbulent field of refractive index, (2) the telescope 
optics, possibly with static aberrations, (3) the AO system optics, (4) the science instru- 
ment optics, that can be merged with the telescope optics. Splitting the phase aberration 
into a static, constant part Tp and a turbulent, zero mean part S(p, 

(p{T,t)=Jp{r) + 6(p{r,t) (1) 

where r is the position vector in the pupil plane and t is the time, and remembering that 
the OTF is also given by the autocorrelation of the phasor exp (iip) in the pupil plane [13] , 
we get, for the long exposure OTF (averaged over an infinite number of realization of the 

2 rcmcmbcr that OTF and PSF are Fourier transforms of each other 
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random AO corrected turbulent phase) 
OTFsys(^) = j- JJ (exp{i[6cp{r,t)-8tp(r + p,t)]}) t x 

P R 2 

exp {i[(p{r) -Tp(r + p)]}P(r)P(r + p) d 2 r (2) 

where v is the angular frequency vector in the focal plane, associated to the spatial shift 
p = Xu in the pupil plane (A is the optical wavelength), S p is the pupil area and P(r) is 
the pupil transmission (1 inside the pupil, outside), and (•) indicates a time or ensemble 
average. Assuming a Gaussian statistics for the phase aberration, which is a very good 
assumption for the uncorrected as well as for the AO corrected phase, it is shown in 
Roddier [H] that the average can be moved into the exponential argument, and we get 

OTFsys W = y J J exp [ - 1/2 D„(r, p)} X 

P R 2 

exp {i[Tp(r) - Tp{r + p)]}P(r)P(r + p) d 2 r (3) 

where 

D„(r, p) = ([Mr, t) - Mr + P, t)] 2 ) t (4) 

defines the phase structure function, a measure of the variance of the phase difference 
between two points separated by a vector p in the pupil. We see that the structure 
function depends not only on the separation vector p, but also on the location r where 
the phase difference is measured. Therefore, if we want to compute the long exposure 
OTF, we need a model in (r,p) of the structure function, which is not necessarily difficult 
to obtain, analytically, but what is more annoying is that we need then to perform a 
numerical integration of Eq. over r for each angular frequency u — p/A. This can 
be a very time consuming effort, and goes against the very objective of the synthetic 
approach. 

Now, it is demonstrated in [H] that the optical turbulence phase - before AO correction 
- is stationary over the pupil, i.e. its statistical properties do not depend on the location 
r inside the pupil. The corrected phase, on the other hand, is not stationary, and its 
residual variance increases from the center to the edge of the pupil. This being said, 
this non stationarity affects mostly the first orders - piston, tip-tilt, defocus ... - and 
the corrected phase can be considered to be stationary for a moderate to high order AO 
system (i.e. moderate to high Strehl). If the phase is stationary, which we will assume 
from now on, the phase structure function can be written D v (r, p) = D^(p) and the 
structure function exponential can be extracted from the integral in Eq. ([3]), so we get 

OTFsysW « exp [ - D„(p)/2] i- J J exp {z[p(r) - Jp{r + Af)] }P(r)P(r + Af) d 2 r 

P R 2 

= OTF AO (^)OTF TS cW (5) 
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The system's OTF can therefore be seen as the telescope OTF filtered by an AO OTF 
filter, exp[— D ip (p)/2\. Separating the OTF this way is equivalent to splitting the overall 
optical system into two independent systems: the optical turbulence plus AO system 
on one side, and the telescope plus instrument optics on the other. The first system is 
therefore not related to the pupil optics in any way, and its description does not include 
the pupil boundaries anymore. For this very reason, synthetic modeling is sometime 
referred to as infinite aperture modeling. 

We discuss now the relationship between the stationary SF and the phase s-PSD. Thanks 
to the stationary assumption, the AO system can be considered as an optical system 
applying a correction on a turbulent phase, regardless of any beam boundaries, all over 
an hypothetical plane perpendicular to the telescope optical axis. Everything looks as if 
the phase was pre-corrected by the AO system before being intercepted by the telescope 
beam. Now, it is shown in [T5] that the stationary structure function is related to the 
spatial correlation of the phase, B^, which is itself equal to the Fourier transform of the 
phase s-PSD (written E v here), and we get 

D» = 2 [B(0) - B(p)] = 2 |y [1 - cos (2ni p)] S^f) d 2 f (6) 

The integral of the phase s-PSD gives the phase variance, and the cosine term is there 
instead of the usual complex exponential we have in a Fourier transform, because as the 
phase s-PSD is even, the sine component of the FT would be naturally null. With this 
last equation, we have completed the link between the phase s-PSD and the long exposure 
PSF. 

To summarize, the procedure to get the long exposure PSF from the phase s-PSD is the 
following: 

1. the phase s-PSD is computed from the analytical expressions given in the next 
sections, for each of the wavefront error components, 

2. the phase SF is computed from the phase s-PSD with Eq. ([6]), using a numerical 
Fourier transform algorithm, 

3. the AO-OTF filter is given by the exponential of the SF, exp (— SF/2), 

4. the telescope OTF is computed analytically or numerically and is filtered by the 
AO-OTF, 

5. the long exposure PSF is obtained by applying a numerical Fourier transform on 
the final OTF. 

The whole procedure therefore consists in the evaluation of a few analytical expressions 
and the computation of two numerical Fourier transforms. 
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3 SPATIAL FREQUENCY POWER SPECTRUM OF THE 
AO CORRECTED PHASE 



3.1 The fundamental equation of adaptive optics 

The starting point for the development of the phase s-PSD is the so-called fundamental 
equation of adaptive optics, which states that at any instant t, the residual wavefront 
error w e is the difference between the incoming atmospheric turbulent wavefront w a and 
the mirror command^ w c - see figure [TJ 



where r is the location in the pupil plane, is the angular separation between the science 
object and the guide star (assumed on-axis without loss of generality), A s is the science 
observation wavelength, and A m is the wavefront sensing wavelength. Note that in our 
initial paper, we used the phase instead of the wavefront in the fundamental equation, 
but we believe now that using the wavefront formulation is more appropriate because it 
is actually the wavefront which is corrected in an AO system. We will therefore develop 
equations for the residual wavefront s-PSD, which will be transformed into the phase s- 
PSD by multiplication with the usual factor (2ir/X s ) 2 . Polychromatic PSF will be modeled 
by computing and averaging the phase s-PSD over the chosen optical bandwidth. 

Including air refractivity The air's refractive index depends (slowly) on the wave- 
length, and measuring the wavefront at a different wavelength than the science observation 
channel introduce a small but noticeable error for systems with a tight wavefront error 
budget. Formally, the air's refractivity (N=n-1) is given by the sum of the refractivity of 
the air's constituents (nitrogen, oxygen, water, carbon dioxide etc.) Practically, though, 
it is shown in [TH] that N can be written as the sum of a continuum and anomalous terms. 
The anomalous terms are associated with the excitation/absorption lines of water vapor 
and carbon dioxide (others constituents have a negligible impact), and because the atmo- 
sphere is naturally opaque at theses wavelengths, the anomalous terms are of no interest 
to us. The origin of the continuum term is actually not different from the anomalous 
terms: it is a sum of the wings of the strong nitrogen/oxygen/ozone excitation lines in 
the ultraviolet which extend to the visible and infrared wavelengths. There are very good 
theoretical/empirical models for this continuum [TT|IT8] in the visible and infrared up to 
10 jum, that are function of the air temperature, pressure, humidity and carbon dioxide 
content. We will not dig here into these models. What is of interest for us is that these 
models can be written as the product of a chromatic term which depends only on the 
wavelength, and a non-chromatic term which depends on the other variables (tempera- 
ture, pressure etc.) Therefore, within the transmission windows of the atmosphere (i.e. 

3 the shape of the mirror is actually set to half of the mirror command, because the OPD on a reflecting 
surface is twice the surface error 
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inside the photometric bands), the wavefront error is simply proportional to the air's 
refractivity, and we can rewrite the fundamental equation of AO as 

w e (r, 0, X s , t) = w a (r, 6, A s , t) - u(\ m , A s ) w c (r, X s , t) (8) 

where we define z/(A m , A s ) = N(X m )/N(X s ) as the dispersion factor. The later formulation 
allows us to develop the model of the mirror command w c for a single wavelength - here 
the science wavelength - and correct for the fact that the wavefront is actually measured 
at another wavelength. The effect of dispersion in discussed with more details in Jolissaint 
and Kendrew [T9j . 

3.2 The residual wavefront error in closed loop mode 

Continuous process assumption AO control is a discretized process: the DM shape 
is updated periodically at a loop period At, with a time delay ti ag following the end of 
the wavefront sensor exposure. In-between these updates, the DM shape is kept constant 
(at least classical systems are working this way). AO control is therefore an integral and 
hold process, in a sense that the updated DM shape is equal to the previous one plus 
a weighted estimate of the wavefront residual. Now, as our approach is stationary in 
nature, no specific instant can be privileged: in other words, the equations we are writing 
need to be applicable at any instant, therefore the discrete integral control equation 
Ck = g x ek + Cfc_i, where Ck represents the DM command at instant tk and g x the 
error signal weighted by the loop gain g, needs to be replaced by the continuous integral 
c(t) = g x J"* At e(t)dt + c(t — At). This approximation is equivalent to assume that the 
closed loop control is applied continuously, as if at any instant t, the DM is updated with 
a command computed from the WFS measurement an instant t — t\ ag earlier. 

The deformable mirror spatial transfer function The command applied to the 
DM is made of the projection of the WFS measurement onto the DM space, whose basis 
is the N-dimensional set of the DM influence functions I i=1 jy. The equivalence of this 
projection in our stationary approach is the convolution of the reconstructed wavefront 
with the DM spatial response, or, in the spatial frequency domain, the multiplication of 
the wavefront Fourier transform with the DM spatial transfer function. The DM spatial 
response is defined by the projection of the Dirac impulse onto the influence function basis 

M=l...N, 

N 

7dm (r) = Pi X i( r ) ( 9 ) 
i=i 
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The coefficients pi are computed from the minimization of the quadratic distance between 
the Dirac impulse and its projection, and we find 
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where S is the covariance matrix of the DM influence functions. 

Note that as the DM response, by definition, is not supposed to vary across the DM 
surface, while actually the projection of the Dirac impulse does depend on its location 
u within the actuators grid, we will simply consider the average DM response over all 
possible locations u. We have found that this ad-hoc procedure generates DM transfer 
functions that better represent transfer functions measured on real systems. Practically, 
as the DM actuator array is periodic, we find that hi as the average of the influence 
function X over the square space □ centered on the optical axis and of width equal to the 
inter-actuator pitch. For another actuator grid geometry, the space over which the DM 
response is averaged would be different. 

The DM spatial transfer function is given by the Fourier transform of the DM spatial 
response, and as the z-th influence function at a position rj can be written as the central 
influence function Xq shifted by — r^, we find 



where Zo(f) is the Fourier transform of the central influence function. There are numerous 
models for the DM influence function, all depending on the type of DM technology - see 
for instance [20J. We have developed our own empirical models for a Xinetics, Inc. 177 
actuators DM model and a Boston MEMS 144 actuators model, with a modeling error of 
less than about 0.1 % in amplitude, and used these to compute the DM transfer function. 
A cut through these two DMs response and transfer functions is shown in Figures [2] and 
is compared with Gaussian and pyramid influence function models. 

It is important to note that any other DM basis of independent functions can be used to 
build the DM transfer function: the choice we made here of using the influence function 
basis was dictated by the fact that we had at our disposal accurate empirical models of 
several influence functions types. Now, as pointed out by a reviewer of this paper, influence 
functions are not always a good choice to model the DM surface, and we certainly agree 
that it is particularly true for low order modes (for instance tilt is badly represented by an 
addition of influence functions). If correction of the low order aberrations is of particular 
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Figure 2: Left: Profile of the DM response for a Xinetics (Inc.) and a Boston MEMS de- 
formable mirror, compared to the response function for a Gaussian and pyramid influence 
function model. Right: DM spatial transfer function power (modulus square) associated 
with the DM responses. 

concern for the modeling of a particular system, it might therefore be worth to use a basis 
particularly adapted to the representation of these modes (for instance, as proposed by 
the reviewer, a DM surface generated by a cubic spline interpolation). 

Now, let us develop the DM command w c , including the continuous assumption and the 
DM response. The DM command w c at any instant t is the update of the previous DM 
command (an instant At earlier) with the residual wavefront error, weighted by the loop 
gain g loop 

w c (r, A s , t) = g loop 7 DM (r) * K{m(r, X s , t) + n(r, t)} + w c (r, X s , t - At) (14) 

where * is the convolution product, 1Z is the operator associated with the wavefront 
reconstruction from the WFS measurement m, and n is the WFS measurement noise. 
Note that this equation is not specific to any type of WFS nor reconstruction operator. 

The WFS output m is defined by the measure (direct, gradient or laplacian) of the wave- 
front residual error w e averaged over the WFS integration time At, delayed in time by the 
lag ti due to the WFS readout time and the command computation. Using the notation 
q to indicate a time average with a time lag, we write the WFS measurement as the 
application of a wavefront measurement operator M. on the instantaneous residual error 
w e in the direction of the NGS, which is assumed on-axis (i.e. 0=0), 

m(r, X s , t) = M{ur e {r, 0, A s , t)} = M{w^(r, 0, A s , t) - uT c {y, X s , t)} (15) 
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Separating the deformable mirror space and the wavefront analysis space For 

further developments, we need now to split the atmospheric wavefront into the component 
which is corrected by the DM, and the component which is simply reflected off the DM 
surface, 

W a = Wg *7dm + Wg * (5 - 7 DM ) (16) 
corrected by the DM reflected by the DM 

Beside, the WFS samples the wavefront with a spatial interval A WPS equal to the lenslet 
separation distance. This naturally split the spatial frequency domain into a low spatial 
frequency domain, i.e. the frequencies that can be seen by the WFS, and therefore 
corrected, below the WFS spatial cutoff frequency 

/ WFS = 1/(2A WFS ) (17) 

and a high spatial frequency domain, above / WFS . For a Shack-Hartmann WFS, the lenslet 
array has a square geometry, therefore the low spatial frequency domain is defined by the 
inequalities \f x \ < f WFS and \f y \ < f WFS , a square, and the high spatial frequency domain 
by the complement of this square. One finally gets, with Eq. ffT6l) . 

W a = (MV LF + MVhf) * 7dm + (»a,LF + 

The later formulation allows the separation of the DM actuators and WFS lenslet array 
grid architectures, which can be now studied independently. 

From this point, our model development is done in the spatial frequency domain only. 
Our final objective is indeed to write an equivalent spatial frequency power spectrum filter 
for the AO correction. In the spatial frequency domain, the fundamental equation of AO 
becomes, where q indicates the Fourier transform of q, 

w e (i , 0, X s , t) = w a (f, 0, X s , t) - v(X m , X s ) w c (f, X s , t) (19) 

the DM command becomes 

w c {f , A s , t) = <7 loop (f) r DM (f)£{m(f , A s , t) + n(f, t)} + w c {f, A s , t - At) (20) 

and the WFS measurement 

m(f , A s , t) = M{w a (i, 0, A s , t) - w c (f, X s , t)} (21) 

Inserting Eq. (j21j) into Eq. fl20l) . with Eq. fTT8|) . and as the reconstruction and mea- 
surement operators become spatial filters in the frequency domain, it comes (we drop 
momentarily the common variables to shorten the notation) 

w c =g loop r DM (l - T DM )KMw a , LF 
+ g loop Tl M nM 
+ fi'ioop r DM (i - T OM )TZM 
+ g loop Tl M TZ M w atHF 
- fi'ioop r DM 7^ M w c 
+ fi'ioop r DM 7?. h + w c (t - At) 
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In our approach, the product r DM (l — r DM ) would describe the projection onto the or- 
thogonal of the DM space, followed by the projection onto the DM space. This product 
naturally has to be replaced by the null operator. Beside, the wavefront reconstruction 
action is to revert the WFS measurement, therefore, for a wavefront g LF strictly limited 
to the WFS low frequency space, the cumulated operation of wavefront measurement and 
wavefront reconstruction is the identity operator, i.e. TZ{A4{q LF }} = q LF . With the later 
remarks, and noting that the DM command, by nature, belongs to the low frequency 
WFS space, Eq. ((22} simplifies to 

W C — S'loop J- DM W Cl,LF 
<?loop J- DM W c 

+ fl-ioop r DM "^n + w c (t - At) 

A note on the loop gain Nothing prevent us from setting the loop gain g loop here as a 
free variable in f , too. This allows optimization of the loop gain frequency-by-frequency, 
which is equivalent - in our stationary approach - to modal gain optimization. We will 
therefore keep the notation (?i oop (f ) even if loop gain optimization is not discussed further 
in this paper. 



Independence of the turbulent layers and Taylor hypothesis In the atmosphere, 
optical turbulence is distributed in thin independent layers, each being characterized by 
(1) the so-called refractive index structure constant Cj fi , a measure of the variance of the 
refractive index spatial fluctuation within the layer, (2) the apparent velocity of the tur- 
bulent layer - see for instance [21] . As seen from the pupil of the telescope, the time scale 
over which the wavefront associated to each layer evolves significantly is generally longer 
than the time it takes for the wind to push the layer across the telescope beam. Therefore, 
in first approximation, everything looks as if the optical turbulence profile was made of 
a certain number of frozen wavefront screens translating across the telescope beam with 
the layers wind speed and directions. This assumption of frozen optical turbulence layers 
is referred to as the Taylor hypothesis in the literature^, and has the nice consequence 
that it is possible to transpose, within the layer, a shift in time At into a shift in space 
Ar = vAt, with v the layer's wind velocity. 



Assuming independence of the turbu 
summed over the iVj turbulent layer: 



. ent layers, the correction of the total wavefront 
is equivalent to the correction of each wavefront 
from each layer taken individually, as in any case, cross terms between the layers will 
vanish on average in the computation of the long exposure residual phase s-PSD. Let us 



4 departure from the Taylor hypothesis is discussed in [33] 

5 about 10 layers are generally needed to model a turbulent profile 
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therefore compute the wavefront error spectrum w e j associated with the layer 

As shown in [6], in the Fourier domain, the time average of a wavefront q over a time 
interval At, followed by a time lag t\ ag becomes, for the layer I with a wind velocity Vj, 

|,(f , t) = sinc(At f ■ Vl ) exp [2m (At/2 + t lag ) f ■ vj] (24) 

where At/2 + ti ag , equaling the time interval between the middle of the WFS exposure 
time and the application of the new DM command, represents the overall time lag. With 
the later, the DM command spectrum Eq. (1231) becomes, for the layer I, 

w C: i = g loop Tl M smc(At f • v z ) exp [2ni (At/2 + t iag ) f ■ v,] u7 a , LFi /+ 

5f loop r^ M sinc(At f • v;) exp [2m (At/2 + t lag ) f • vj] KM w a ,HF,i- 

_ \2o ) 

fi f ioo P r DM sinc( At f ■ vj) exp [2m' (At/2 + t lag ) f ■ v z ] w c j+ 

exp(2ni Atf ■ vi)w c ,i 

where we have replaced the time shift At within w c j(t — At) by its equivalent phase change 
exp (2iri Atf ■ vi) in the Fourier domain, using the Taylor hypothesis. 

The noise term, being added to the wavefront slope measurement, is not linked in any way 
to the turbulent layers, therefore it does not make any sense to write a noise term for each 
layer. The noise term is independent from the other error terms and needs to be treated 
separately, so we do not include any noise term in the last equation. The overall WFE 
will simply given by the sum of the servo-lag contribution and the noise contribution. 

Regrouping the terms in w Cs i, we end up with an expression for the DM command 

w Cil (f , X s , t) = [g loop (f) T 2 DM (f) sinc(At f ■ v t ) exp [2m (At/2 + t lag ) f • v,] w a , LF ,,(f , 0, A s , t)+ 

g loop (f) r| M (f ) sinc( At f ■ v,) exp [2m (At/2 + t lag ) f • v,] £(f ) M(f) w ml (f, 0, A S) t) }/ 

1 + g loop (f) T DM (f) sinc(At f • vj) exp [2m (At/2 + t lag ) f ■ vj - exp (2m At f ■ v f ) 

(26) 

Inserting Eq. (1261) into Eq. (1191) . we get, again for the layer I, 

W e ,l(f, 0, X s , t) = U7 a ,HF,i(f , 0> K, t) + U> a ,LF,i(f, 0, X s , t) 

v(X m , A 8 )ff loop (f ) T^ M (f ) sinc(At f ■ vQ exp [2m (At/2 + t lag ) f ■ Vj]w a , LF , t (f , 0, X s , t) 
1 + # loop (f ) T DM (f ) sine (At f ■ v,) exp [2m (At/ 2 + t lag ) f ■ v,] - exp (2m At f ■ vj) 

^(A m , A a )a„ op (f ) ryf) sinc( At f ■ v,) exp [2m (At/2 + t lag ) f ■ vjggMg^ 0, X 8 , t) 
1 + # loop (f ) T DM (f ) sinc(At f ■ vi) exp [2vrz (At/2 + t lag ) f ■ vj] - exp (2™ At f • v,) 

(27) 

As an angular shift 6 is seen, at the layer altitude hi, as a spatial shift Ar = h\Q, using 
the shift theorem of the Fourier transform we get 

«VF,z(f , 0, A s , t) = exp (2m h t f ■ 0) w a ^ ¥) i(f, 0, A s , t) (28) 
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We can now rewrite Eq. (J2j 

w e> i(i, 0, X s ,t) = 



V 

high order WFS "fitting" error 



(f,0,A„t) (29) 



aniso- servo error 



+ ^AL,l(f) (f,0,A a ,t) 

V v ' 

WFS aliasing error 

and identify the four fundamental terms of the residual wavefront error: 

1. the high order WFS error, usually named the "DM fitting error" in the AO literature 
- which is actually for us the part of the atmospheric wavefront which is not seen 
by the WFS, therefore cannot be corrected by the DM, so we think that calling this 
error the high order WFS error is more appropriate, 

2. the angular anisoplanatic AND loop servo-lag error, identified here as the "aniso- 
servo" error, as anisoplanatism and servo-lag error are correlated (with the Taylor 
hypothesis, a wavefront shift in time can be compensated by a negative wavefront 
shift is space), 

3. the WFS aliasing error: the high order wavefront error is seen by the WFS as a 
low spatial frequency error and reconstructed as such, therefore the AO system is 
compensating an error which actually is non-existant, 

4. and finally the WFS noise term, discussed later. 

F AS i is defined as the aniso-servo spatial filter for the layer I, 



(30) 



F AStl (f) =exp (2m fyf -0)- 

A 8 )fl loop (f ) Tg, M (f) sinc(At f ■ vQ exp [2m (At/2 + t lag ) f ■ vj 
1 + g loop (f) T DM (f ) sinc(At f • vj) exp [2m (At/2 + t lag ) f • Vl ] - exp (2m At f • v,) 

and F AL j is the WFS aliasing spatial filter for the layer I 

KA m , A 8 )fl loop (f ) Tg, M (f ) sinc(At f ■ v;) exp [2m (At/2 + t lag ) f ■ vjg(f)jf(f) 
" 1 + g loop (f) T DM (f ) sinc(At f • v,) exp [2m (At/2 + t lag ) f ■ vj - exp (2m At f • v,; 

It is interesting to examine the limits of the aniso-servo and WFS aliasing spatial filters 
when there is no angular separation between the science object and the NGS, i.e. 6=0, 
and when the WFS integration time and loop lag are set to zero, Ai=ti ag =0. We find 

f A s,i — > i - v r D M 

F^, -)• -uT DM nM 



(31) 
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which indicates that in the absence of aniso-servo error, the residual low frequency wave- 
front error is generated by (1) the refractive index dispersion - and we see as we would 
expect that this error is proportional to the dispersion factor, and (2) the aberrations seen 
by the WFS that the DM cannot correct, which are never null because even if the DM 
actuator pitch is equal to the WFS lenslet pitch, a perfect correction of the low spatial 
frequencies would require a DM with sinus cardinal influence function (Fourier transform 
of a sinus cardinal is a door function), which is only approximated by actual influence 
functions - see Figure [2J The same is true for the WFS aliasing error. This behavior 
corresponds well to what we would have expected. 

A note on the correlation between the error terms Computing the s-PSD as- 
sociated with the four fundamental wavefront error terms above essentially consists in 
computing the modulus square of Eq. (|2T?|) . averaged over the time t. Therefore, cross 
products appear between the four error terms. Now, the noise term is naturally not cor- 
related with the other errors and can be treated separately. The high order WFS error 
is not seen by the system and transmitted to the output of the system, unaffected. We 
will assume in this paper, without discussing it further, that the cross products between 
the low and high order spatial frequencies are negligible relative to the main error terms. 
Therefore, in what follows, the four terms above will be discussed independently from 
each others. 

We will now make use of the expressions developed above for the wavefront error terms 
to develop the analytical expressions of the residual phase s-PSD of the four fundamental 
errors we have identified. As the residual phase s-PSD is given by spatial filtering of the 
optical turbulence phase s-PSD, we will start by recalling the expression of the later, as 
given in the literature. 

3.3 The turbulent phase spatial power spectrum 

The s-PSD of the turbulent phase is discussed in Roddier p3]. In the atmosphere, the 
extension of optical turbulence is necessarily limited by the individual layers thickness, 
and an optical turbulence outer scale L was included in the Kolmogorov s-PSD to account 
for this spatial limitation. This modified Kolmogorov s-PSD is called the von Karman 
s-PSD in the literature (see for instance Winker [23], and Maire et al. [23] for a few other 
s-PSD models) and is given, at the wavelength A, by 

~ ATM (f) = 0.0229 r (A)~ 5/3 (|f | 2 + 1/^)~ U/6 (32) 

where vq is the Fried parameter, a measure of the strength of turbulence, defined as the 
telescope diameter whose focal plane angular frequency cutoff would be the same than 
the optical turbulence cutoff frequency (see Fried [23]). r is generally given at 500 nm 
in the literature, and we will follow this convention, unless indicated differently. Typical 
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values for r$ at 500 nm extend from 5 cm (bad observation site, day-light conditions) 
to 25 cm (excellent site). The optical outer scale is generally in the range 20 to 40 m, 
surprisingly with very few variations between the different sites where this quantity has 
been measured. 

The Fried parameter is associated to the vertical profile of the optical turbulence structure 
constant C%(h), following 



r 



POO 

(\y 5/3 = 0.4234 (27r/A) 2 / C 2 N (h)dh (33) 

Jo 



which can be written, in the case of Ni independent optical turbulence layers, as a sum 

r (A)" 5 / 3 = 0.4234 (2tt/A) 2 £ C 2 N>1 Ah t = >T r ,(A)^ 5 / 3 (34) 

i=i i=i 

where Ahi is the layer's thickness, and r 0j / defines the layer's Fried parameter. Conse- 
quently we can also define a phase s-PSD for each layer, 

"atm, 2 (f ) = 0.0229 r XA)" 5/3 (|f| 2 + 1/^)~ 11/6 (35) 
which naturally sums up to the overall phase s-PSD, 

Ni 

"ATM (f ) = H ATM) j(f) (36) 

1=1 

Using this notation, the phase s-PSD associated with the low and high order turbulent 
wavefront errors - io a ,LF,i and w a ^ HF j - will be written now 

(27r/A s ) 2 (K LFii (f , 0, A s , t) \ 2 ) t = /i LF (f) S ATM ,,(f) (37) 
(27r/A s ) 2 (|^ aiHF Xf,0,A s ,t)| 2 ) i = y u HF (f)H ATMi/ (f) (38) 

where /i LF and low and high spatial frequency masks, defined, for a Shack- 

Hartmann WFS, by the square domain 

/i LF (f) = < (39) 
I elsewhere 

and ^H*(f) = l-M f ) ( 4 °) 

3.4 The high order WFS spatial power spectrum - or "fitting error" 

The high order WFS phase error s-PSD is simply given by the atmospheric turbulence 
phase s-PSD limited to the high spatial frequency domain and is written 

~ HF (f) =/i HF (f)S ATM (f) (41) 

where S ATM is given by Eq. 
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3.5 The aniso-servo spatial power spectrum 

The aniso-servo phase error s-PSD is given by the time average of the modulus square 
of the aniso-servo wavefront error, translated into a phase error, summed over the Ni 
turbulent layers. From Eq. f[2"§j) and Eq. f[3"T|) . we get 

2 AS (f) = // LF (f) l F As,d 2 (f) 2 ATM ,,(f) (42) 

where S ATMj ; is given in Eq. (135]) and 

|F AS /(f) = (l+^ p (f)rL(f)sinc 2 (Atf-vO[l + ^(A m ,A s )rL(f)]/2-cos(2 7 rAtf-v i ) 

+ # loop (f)r 2 M (f) sinc(Atf • vi)u(X m , A s )x 
{ cos [2-irhif ■ 6 + 2vr(At/2 - t lag )f • vj] - cos [27i7ijf ■ - 2vr(At/2 + t iag )f ■ v,]} 
+ ^oo P (f)r DM (f) sinc(At f • v,){ cos [2vr(At/2 + t lag )f • v,] - cos [2vr(At/2 - t lag )f ■ v ; ]} 

- ^oo P (f)rL(f) sinc 2 (Atf ■ v,MAm, A.) cos (vr/^f • 0)) / 

(l + 5l 2 oop (f)r 2 M (f) sinc 2 (At f ■ v,)/2 + <? loop (f)r DM (f) sinc(Atf • v,)x 

{ cos [2tt( At/2 + t lag )f • v,] - cos [2tt( At/2 - t lag )f • vj] } - cos (2?rAtf ■ v z ) ) (43) 

While this equation seems impressive, coding it into a computer program does not repre- 
sent a particular challenge. 

3.6 The WFS aliasing spatial power spectrum 

The aliasing error is given, for the layer I, by the term - see Eq. (T2~T|) 

wam(M) = -n(f)M(f)w a>HF>l (f,o,x s1 t)x 

KA m , A B )fl loop (f ) T 2 M (f ) sinc(At f ■ v;) exp [2ttz (At/2 + t lag ) f ■ vj 
1 + & oop (f ) T DM (f ) sinc(At f • v,) exp [2m (At/2 + t lag ) f • vj - exp (2m At f ■ vj) 1 J 

the product 7?. .M uj ,hf,z was already developed in our initial paper, and is recalled here. A 
Shack-Hartmann WFS produces a measurement of the wavefront slope (gradient) in both 
x and y directions, with a spatial sampling given by the WFS lenslet array pitch A WPS - 
which is also the lenslet width. In the spatial frequency domain, the slope measurement is 
given by the multiplication of the wavefront Fourier transform with the two components 
operator (one for each direction) 



M(f) = [M x (f),M y (f)} = 27rzA 2 VFS [fsinc(A WFS / x )sinc(A WFS / J/ )] * III(A WFS f) (45) 

where the product with the spatial frequency vector f stands for the derivative in the 
Fourier domain, the sine function is for the wavefront average over the lenslet area, and 
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III(A WFS f), the Dirac comb, represents the recurrence of the measured slope spectrum 
with a spacing 2/ WFS in both x and y directions, and is responsible for the aliasing of 
the part of the spectrum above the WFS cutoff frequency / WFS inside the low spatial 
frequency domain - which always occurs, because the turbulent wavefront spectrum is 
not band limited at high spatial frequency. 

The analytical expression for the reconstruction operator Fourier transform 1Z is computed 
from the minimization of the quadratic distance between the slope measurement and the 
actual slope (least square error algorithm, LSE). The weakness of the LSE algorithm is 
that the WFS noise is reconstructed as a real signal, without penalty. Other algorithms 
have therefore been proposed that make use of a priori knowledge of the Kolmogorov- 
statistics based signal and noise statistics to minimize the contribution of the noise on 
the reconstructed signal - see for instance [12]. A discussion of the pros and cons of these 
different algorithms is beyond the scope of this paper, though, so we will stick with the 
LSE-based algorithm, as it is the most simple and most straightforward to implement. 

As we have seen, the slope measurement operator in the Fourier domain is, basically, a 
multiplication with the spatial frequency vector f . The reconstruction operator is therefore 
the inverse operator, i.e. the inverse of the vectojfjf, but ignoring the Dirac comb, because 
the reconstruction does not extend beyond the WFS cutoff frequency / W fs> and we find 
(the factor A WFS disappears because it is actually part of the Dirac comb convolution 
product) 

m - fern, *.(*)] - 2xi | f |, sinc(A J / j siD c(A w , s /,) < 46 > 

We can now develop the term TZM w a ,HF,i from Eq. (jUJ). Using the two equalities 

HI(A WFS f) = — ^-^-Jv-l^-) (47) 



and 



we find 



7V WFS m ^ A WFS A WFS 



sinc(AwFs/[a;,y] - [m,n]) _ (-l) [m ' n] A WFS f [Xiy] 
smc(A WFS f [Xiy] ) A WFS /[ Xiy ] - [m,n] 



(4£ 



ft(f)A*(f)< HF ,,(f,0,A S) t) = 

E (-!) m+n (tt^ + tt^W/-- ir->fy- ^->M.,t) (49) 

m,n=-oo \JV Awfs Jx A W FS / AwFS AwFS 

\m\ + \n\>0 

where it is important to note that \m\ + \n\ > because we do not want to include the 
low frequency part of the wavefront spectrum in the sum (as it is of course not aliased). 



' • v = 1 has the solution u = v/\v\ 2 
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We can give now the expression for the WFS aliasing phase error s-PSD. Summed over 
the Ni independent layers, it is given by 

oo 

S AL (f) = E "AL,/(f) (50) 

where, for each layer, 

S AM (f) = (|^ AM (f,t)| 2 )* = ^(A m ,A s )^ oop (f)r 4 DM (f)sinc 2 (Atf • v,)/ 

- ^ p(f)rL(f) sinc 2 (At f • v,)/2 + <? loop (f)r DM (f) sinc(At f ■ v,)x 
g )f ■ Vl ] - cos [2vr(At/2 - t lag 

m „ n 



{cos [27r(At/2 + t lag )f ■ vj] - cos [27r(At/2 - t lag )f ■ v { ]} - cos (2?rAtf • x 

lf|4 f n + r m I "ATM,IUa: a j a J l di J 

I 1 ! m,n=-oo \ /y Awfs Jx A W fs / AwFS AwFS 

|m|+|n|>0 

where we made the assumption that the correlation of the phase for frequencies sep- 
arated by a 2/ WFS interval is negligible. It is important to realize that the term 
T = ^j?E(-) S ATM,i(-) m Eq. (P) has singularities at (f x = 0,m = 0), (/„ = 0,n = 0) 
and f=0. Computation of the limits gives 



^(/an 0) = ^m=-oo ^ATM,l(fx — A ™ j 0) 
|m|>0 WFS 



(52) 



3.7 The WFS noise spatial power spectrum 

From Eq. 02 3p . it comes 

H NS (f) = (|^ NS (f,t)| 2 ) = ^ 2 (A m ,A s )^ oop (f)r 2 M (f) (|£(f)n(f,t)| 2 > (53) 

What is the noise n(r,t) made of, exactly ? it is a discrete quantity, in space and 
time, made of two components in x and y, n(r,t) = [n x (r , t) , n y (r , t)] sampled on a 
grid of spacing A WFS . Its spatial spectrum is therefore necessarily limited to the do- 
main \f x \, \f y \ < 1/(2A WFS ). As the noise over the lenslets is uncorrelated, all values are 
possible at any instant and location, therefore the noise s-PSD is necessarily white, i.e. 
it is constant in the domain \f x \, \f y \ < 1/(2A WFS ), and is the same for both x and y 
components. So, we can define 

(\n x {f,t)\ 2 ) t = (\h y (i,t)\ 2 ) t = AT 2 (f) = constant (54) 

such that 

A7 2 (f)d 2 / = ( x 2 EA . CL =A7 2 /A 2 VFS (55) 



l/-l.l/»l< 

1/(2A W fs) 
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A/" 2 =AwFS^EA,CL (56) 

where <r 2 EA CL is the closed loop noise equivalent angle (NEA) variance, discussed later. 
The noise s-PSD is given by the average modulus square of w NS . With the reconstructor 
Fourier transform - given in Eq. ( |46|) . it comes 



(|^(f)n(f,t)| 2 ) = <|^(f)n,(f,t)| 2 > + (\TZ y (f)n y (i,t)\ 2 ) 



f 2 A\n^M) + f 2 y(\n y (i,m A^a 2 



Air 2 |f | 4 sine (A WFS f x ) sine (A WFS f y ) Air 2 |f | 2 sine (A WFS f x ) sine (A WFS f y ) 
therefore we get, from Eq. ( 15T?j) . 



(57) 



— ,n\ V (^rni A s ) T DM (f ) A wpg <7 NEAi CL 

" NSI ' ~ 4vr2|f|2 si nc 2 (A WFS / x )sinc 2 (A WFS /,) 1 } 

Note that the loop gain does not appear anymore directly in the noise s-PSD. Indeed, the 
noise s-PSD has to be seen as a spatial filter, actually not different from its formulation in 
open loop, but where the NEA is now a closed loop NEA. In other words, it is the NEA 
which is affected by the closed loop noise transfer function, not the spatial properties of 
the noise s-PSD. Note also that there is no analogy between the spatial frequency white 
noise and this open loop white noise. Indeed, any type of temporal spectrum is possible 
for the NEA signal on the lenslets, and as the lenslet noise is decorrelated from a lenslet to 
another, all noise distribution have the same probability over the lenslet array, therefore 
the spatial noise distribution is white whatever the lenslet noise statistics. In other words, 
it is the independence of the noise from a lenslet to another which enable the separation 
of the noise s-PSD from the noise temporal PSD. 

We discuss now the closed loop NEA variance. Let us consider a classical model of the 
loop architecture: a WFS with an integration time At, followed by a delay t\ ag due to 
the WFS readout time and the command computation time, then a numerical integral 
controller with gain g loop , a digital to analog converter, and the DM. The noise rejection 
temporal transfer function, defined as the ratio between the NEA signal and the residual 
error, is given, in the steady state, by (see Demerle et al. [26] for a detailed discussion), 

H , s = 0ioop(f ) exp (-27riti ag u) 2iriAt v 

n[U) -(27rAtz/) 2 + ^ loop (f)exp(-27rit lag z/)[l-exp(-27T2Atz/)] 1 ' 

where v is the temporal frequency. The noise power transfer function is given by the 
modulus square of H n , and we find 



H n (v) 



g loop (i) 2 (2ttA*i/) 



(27rAtz/) 4 -2(7 loop (f)(27rAtz/) 2 {cos(27rt lag z/)-cos [27r(At+t lag ) i/]}+2^[i_cos(27rAtz/)] 

(60) 
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The closed loop NEA variance cr 2 EA CL is given by the integral of the filtered open loop 
temporal power spectrum, and as the later is a white noise limited to the domain \u\ < 
1/(2 At), we find 

l/(2At) 

a 2 NEAiCL = 2Ata 2 NEAiOL J \H n {v)\ 2 dv (61) 

o 

The open loop NEA variance cx^ EA OL depends on the number of NGS photons received 
per lenslets during the WFS integration time, the WFS geometry (lenslet width), the 
WFS integration time, the WFS detector read noise, and the NGS image size, which is 
tilt-compensated for a closed loop system. Several models have been developed in the 
literature for this term, and will not be reproduced here - see for instance Rousset [27] 
and Thomas et al. [28]. 

Our closed loop phase s-PSD model is now complete, and is given by the sum of the s-PSD 
of the four fundamental AO errors: the high frequency WFS error, the aniso-servo error, 
the WFS aliasing error, and the WFS noise error, 

H v (f) = ~ HF (f) + S AS (f) + E AL (f) + H NS (f) (62) 

Let us now illustrate the usefulness and usage of the synthetic model with a few examples. 
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Table 1: Optical turbulence profile used in our illustrative examples. Paranal observatory 
type (taken from "E-ELT AO design inputs: relevant atmospheric parameters", ESO 
document E-SPE-ESO-276-0206, except for the wind direction, which is set arbitrarily) 



height 


CfjAh 


wind 


wind 


above 


distr. 


speed 


dir. 


pupil m 


% 


ms 


/ x-axis 


42 


53.28 


15 


38° 


140 


1.45 


13 


34° 


281 


3.5 


13 


54° 


562 


9.57 


9 


42° 


1125 


10.83 


9 


57° 


2250 


4.37 


15 


48° 


4500 


6.58 


25 


-102° 


9000 


3.71 


40 


-83° 


18000 


6.71 


21 


-77° 



Table 2: Synthetic model parameters values used in our illustrative examples. Optical 
turbulence parameters are given at 500 nm. 



telescope diameter 
seeing angle 
Fried parameter 
outer scale 
phase time scale 
isoplanatic angle 
dispersion factor 
DM conjugation 
DM pitch 
DM actuator geometry 
DM influence function 



D=8 m 
w = 0.65" 
ro = 15.5 cm 
L = 25 m 
r = 3 ms 
6 = 2.4" 
v = 0.99 
to pupil 
- free - 
square 

Xinetics, Inc. 



WFS type 
WFS lenslet width 
WFS throughput 
WFS detector noise 
WFS integration time 
loop time lag 
loop gain 
NGS location 
NGS magnitude 
NGS BB temp 



SH 

- free - 
31% 

2 e/px 

- free - 
0.8 ms 

- free - 

- free - 

- free - 
5700 K 



Table 3: wavefront errors RMS for the example discussed in section I4TT1 



high frequency 


47 nm 


WFS aliasing 


25 nm 


aniso-servo 


96 nm 


WFS noise 


142 nm 


total error 


179 nm 
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4 ILLUSTRATIVE EXAMPLES 

The synthetic method has been coded into our AO modeling code PAOLA0, a general 
purpose IDL-based toolbox for modeling the AO correction of segmented telescope static 
and optical turbulence aberrations. It includes open and closed loop single NGS mode, 
and a complete multiple NGS ground layer AO mode. We present in this section several 
studies undertaken with PAOLA, to illustrate the usefulness and usage of such a synthetic 
tool. An optical turbulence profile and standard telescope and AO system parameters are 
defined in the Tables [T] and H] to be used in the different examples. 

4.1 Structure of the residual phase spatial power spectrum, and its impact 
on the PSF wings 

We consider in this example a DM with a square actuator grid, an actuator pitch of 20 
cm (as projected in the telescope primary mirror), a SH-WFS lenslet width of 20 cm (as 
projected in Ml), a WFS integration time of 2 ms, a loop gain of 0.5, an off-axis NGS 
at 3" and a NGS magnitude mV=12 (other parameters are given in the Tables [1] and |2]) . 
The wavefront RMS error for the four classical components for this example are given in 
Table El 




pupil plane spatial frequency [1/m] focal plane x — cordinate [asec] 

Figure 3: Profile of the four fundamental wavefront errors s-PSD, for the case discussed 
in section I4.ll HF is for high frequency error, AL for WFS aliasing, AS for aniso-servo 
and NS is for WFS noise. The vertical dotted lines show the transition low/high spatial 
frequencies, at 2.5 m -1 . Right: Profile of the PSF associated with the s-PSD, at 1.25 /jm. 
The spatial frequency and the angular coordinate are at the same scale in both figures, 
i.e. x = \f. The transition core to halo occurs at an AO radius of 0.64" . 
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The s-PSD profile and the corresponding PSF profile are shown in Figure [3j It is well 
known [291130] , as can be seen with this example, that the PSF wings structure mimics the 
PSD shape. Figure H] shows the four fundamental errors s-PSD in the spatial frequency 
plane. The noise and aniso-servo errors affects the central parts of the low frequency 
domain. Anisotropy of the aniso-servo error is a combined consequence of the wind 
direction and the off- axis location of the NGS. WFS aliasing, as it is expected, affects the 
highest spatial frequencies of the low frequency domain. 




Figure 4: s-PSD of the four fundamental wavefront errors, for the case discussed in section 
14.11 Top-left: high frequency error. The central black square shows the low spatial 
frequency domain inside ±/ WFS (here equal to ±2.5 m _1 ). Top-right: WFS aliasing error. 
Bottom-left: aniso-servo error; the s-PSD is elongated in a direction which is a composition 
of the main wind direction (+54° relative to the horizontal axis) and the NGS orientation 
(along the x-axis). Bottom-right: WFS noise error. Please note that the low frequency 
s-PSD figures and the high frequency error figure have different spatial scales: the width 
of the low frequency images is 5 m -1 while the width of the high frequency s-PSD image 
is 20 m -1 . 
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One can therefore expect, in general, the noise errors to contribute essentially to a widen- 
ing of the PSF core, the aniso-servo error to affect the PSF wings in the region between 
the core and the transition to the residual seeing halo (due to the high frequency error), 
while aliasing would affect mostly the transition region. In this example, noise clearly 
dominates the PSF structure, though. 

A note on the spatial frequency pixel size We have seen that both aniso-servo and 
aliasing s-PSD equations include cosine functions of products in f • v and f • 0. These 
cosine terms need to be well sampled in the spatial frequency domain when building the 
numerical matrices f x and f y : an under-sampling would lead to an underestimate of the 
wavefront error variance, as the later is estimated from numerical integration of the s- 
PSD, and an incorrect representation of the PSF wings structures. The consequence of 
such an under-sampling is an over-optimist estimate of the Strehl ratio for large off-axis 
NGS angle (the Strehl would saturate above a certain value while it should absolutely 
converge to zero for larger and larger off-axis angles). The same is true for the servo-lag 
error, where the performance would be over-estimated for large WFS integration time 
and/or high wind speed. Practically, our experience with PAOLA shows that the cosine 
terms should be sampled with at least 10 samples over one period. Nyquist sampling is 
by far not sufficient here. 

4.2 Strehl ratio of the four fundamental wavefront errors 

In this section we simply illustrate how the Strehl associated with the four fundamental 
errors varies with the main system parameter associated to each error: the WFS lenslet 
pitch for the WFS high frequency and WFS aliasing error (Figure [5]), the NGS off-axis 
angle for the anisoplanatic error (Figure El left), and the NGS magnitude for the WFS 
noise error (Figure El right). It is worth noting that these curves were built in only a couple 
of seconds of CPU time (iMac computer, 3.06 GHz Intel Core 2 Duo processor). We tested 
also the Marechal approximation, stating that for low to moderate phase variance a't the 
Strehl ratio is given by S = exp(— a^). See the dashed curves in Figures and El We 
find that this approximation is actually excellent for the high frequency and WFS aliasing 
error, relatively good for the WFS noise error, and more questionable for the anisoplanatic 
error - below a Strehl ratio of about 40% in the example given here. 

4.3 Open loop versus closed loop performance 

We claimed in the introduction that open and closed loop systems behave differently in 
dim NGS conditions, and limiting magnitude might be quite different for both modes. 
This claim came from the realization that the noise transfer function are quite different 
in the two cases, as well as the servo-lag error transfer functions. In order to illustrate 
this, we computed the Strehl at 1.25 microns for our standard conditions, and a NGS 
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Figure 5: Left: High order WFS error and WFS aliasing phase variance as a function of 
the WFS lenslet pitch. It is known - as can be seen here - that for a SH-WFS the aliasing 
variance is about l/3rd of the high frequency error. Right: Strehl ratio and WFS lenslet 
pitch. See section I4.2I Dashed curves shows the Strehl computed from the Marechal 
approximation. 
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Figure 6: Left: Angular anisoplanatism Strehl and NGS off- axis angle. Dashed curves 
shows the Strehl computed from the Marechal approximation. Right: Strehl and NGS 
magnitude. Dashed curves shows the Strehl computed from the Marechal approximation. 
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magnitude in the range 4 to 16. 

Initially, we set the WFS integration time fixed at 1 ms, for both modes, and a closed loop 
gain of 0.5. See Figure We find that the servo-lag error is higher for the closed loop 
mode, because the rejection transfer function has a lower bandwidth in closed loop than 
in open loop. The noise error on the other hand is higher in open loop, and this is because 
the noise is basically unfiltered in open loop, while the noise transfer function is a low 
pass in closed loop, filtering the high temporal frequency of the white noise spectrum. For 
given wind conditions, an open loop system can be run faster than a closed loop system 
because the servo-lag error is intrinsically lower. Therefore, we might think that a dimer 
NGS could be used in open loop. This is actually true only for bright NGS, where the 
increased open loop noise error is still low with respect to the servo-lag error. One can 
see for instance in Figure [TJ right, that for a Strehl specification of 0.9 (as it would be 
for an Extreme AO system), the open loop limiting magnitude would be 4 magnitudes 
higher than for the closed loop system. For dimer NGS, though, the increase of open loop 
noise overcomes the decrease of servo-loop error, and the open loop system performs less 
than its equivalent (same WFS integration time) closed loop system. To summarize, for 
a given WFS integration time, open loop systems outperform closed loop systems only in 
bright NGS conditions. 
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Figure 7: Left: Servo-lag and WFS noise error RMS in open and closed loop mode, for an 
exposure time of 1 ms and a loop gain of 0.5. Right: Strehl associated with the previous 
WFE, at 1.25 /iin (OL: open loop, CL: closed loop). 
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As a final experiment, we optimized the WFS integration time and the closed loop gain, for 
each value of the NGS magnitude. See Figure [HJ Optimization has several consequences. 
First, the limiting magnitude gain is very significant, more than two magnitudes in this 
example. Second, it makes the open and closed loop servo and noise error converge: this 
is explained by the fact that the structure of the servo-lag and noise s-PSD are the same 
in both modes (see Figure H]), therefore optimization converge towards the same solution. 
Note that the open loop advantage for bright stars disappears with optimization: the 
closed loop mode performs the same as the open loop mode at any magnitude. Therefore, 
contrary to what the intuition would tell, from a control efficiency point-of-view, we assert 
that there is no advantage of using an open loop rather than a closed loop scheme in NGS 
AO mode. 
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Figure 8: Left: Servo-lag and WFS noise error RMS in open and closed loop mode, same 
conditions than Figure but with optimization of the WFS integration time and loop 
gain. The optimized open and closed loop modes WFE are now basically indistinguish- 
able. Right: Strehl associated with the previous WFE, at 1.25 //m; dotted line: before 
optimization; continuous lines: after optimization. 



5 CONCLUSION 

This paper presents a synthetic modeling method for closed loop astronomical adaptive 
optics, and complements earlier work on open loop modeling using the same approach. 
The concept of the synthetic method and its complementarity with the more classical end- 
to-end modeling approach is discussed extensively. The main advantages of the method 
is that it allows a rapid and direct modeling of the long exposure PSF without going 
through a long and cumbersome Monte-Carlo process Then, we give the detailed analytical 
calculation of the spatial power spectrum of the residual AO corrected phase, as well as 
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the steps to go from the power spectrum to the long exposure PSF, allowing the reader to 
write his/her own modeling code. Dispersion of the air refractive index is included in the 
model, as well as the deformable mirror spatial transfer function. This method has been 
implemented into our AO modeling toolbox PAOLA, and used to study a few illustrative 
examples of the usage of the synthetic method to explore the performance of closed loop 
AO systems. It is found for instance that when optimizing the WFS integration time, 
open loop and closed loop system have basically the same performance (same limiting 
magnitude). Finally, it is important to recall that the foundations of the method do not 
depend on the type of WFS neither on the type of wavefront reconstruction method, or 
control algorithm. Also, the method is in principle not limited to single NGS case but 
can be extended, as it has been done by others, to multi-NGS and multi-DM modes. In 
this paper, though, we have simply considered the case of an AO system with a single 
NGS, for a Shack-Hartmann type WFS, a classical LSE wavefront reconstruction and a 
simple integrator control. 
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